Física Computacional: Tarea 1¶
Profesor: Pablo Benítez Llambay
Integrantes:
- Martín Raguileo Reyes
- Fernando Zamora Carrasco
IMPORTANTE EJECUTAR ESTA CELDA PARA QUE SE EJECUTE EL CÓDIGO SIN PROBLEMAS¶
import importlib
import subprocess
import sys
def install_if_missing(package, name=None):
name = name or package
try:
importlib.import_module(package)
print(f"'{package}' ya está instalado.")
except ImportError:
print(f"Instalando '{package}'...")
subprocess.check_call([sys.executable, "-m", "pip", "install", name])
# Lista de paquetes a verificar/instalar
packages = {
"numpy": "numpy",
"scipy": "scipy",
"matplotlib": "matplotlib",
"IPython": "ipython"
}
for pkg, name in packages.items():
install_if_missing(pkg, name)
'numpy' ya está instalado. 'scipy' ya está instalado. 'matplotlib' ya está instalado. 'IPython' ya está instalado.
import os
import shutil
import numpy as np
from math import sqrt
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation
from typing import Callable, List, Dict
from IPython.display import Video, display
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
Introducción¶
El estudio de los sistemas presa/depredador constituye un área fundamental dentro de la dinámica de poblaciones, ya que permite comprender como interactuan dos especies que dependen entre sí para su supervivencia. El modelo Lotka-Volterra es uno de los enfoques matemáticos más conocidos para describir estas interacciones a través de un sistema de ecuaciones diferenciales no lineales.
Para este proyecto se utilizan distintos métodos numéricos para la evaluación de aspectos clave dentro del modelo clásico Lotka–Volterra, como la influencia en el paso del tiempo, el análisis de la estabilidad de las soluciones, el análisis del error entre los distintos métodos implementados y la identificación de los puntos fijos del sistema. Asimismo, se realiza una exploración de la variación de parámetros para comprender cómo afectan el comportamiento de las poblaciones, se construyen diagramas de fase que permiten visualizar ciclos poblacionales y se discuten posibles extensiones del modelo.
Marco teórico¶
Modelo de Lotka-Volterra¶
Para obtener las ecuaciones que describen la tasa de crecimiento de presas y depredadores, es necesario partir de un análisis aislado de cada una de las poblaciones y posteriormente realizar un analisis conjunto resultado de la interacción entre ambas poblaciones.
Por un lado, en el caso de que no exista una capacidad de carga $K$ en el ecosistema, se presenta un crecimiento exponencial de las presas, determinado por una tasa de crecimiento en ausencia de depredadores $(\alpha)$, descrito por la siguiente ecuación:
$$x(t) = x_{0} e^{\alpha t} = x$$
por lo que su expresión diferencial, es decir, la tasa de crecimiento queda tal que:
$$ \frac{dx}{dt} = \alpha x_{0}\, e^{\alpha t} \;\;\;\Rightarrow\;\;\; \frac{dx}{dt} = \alpha x $$
Ahora, considerando un análisis conjunto, la población de presas se ve afectada por la depredación, lo que provoca una disminución en su número. Este efecto relaciona las poblaciones de presas y depredadores mediante una tasa de mortalidad de las presas debida a los encuentros con depredadores $(\beta)$. Finalmente, la tasa de crecimiento de las presas se expresa como:
$$\frac{dx}{dt} = \alpha x - \beta xy$$
Para la tasa de crecimiento de los depredadores, la población de depredadores ya está limitada por la disponibilidad de presas en el sistema, por lo que, considerando una caso aislado en la que no existen presas, la población presenta una disminución exponencial, determinada por una tasa de mortalidad natural de depredadores por ausencia de comida $(\gamma)$, descrito por la siguiente ecuación:
$$y(t) = y_{0} e^{-\gamma t} = y$$
por lo que sus expresión diferencial, es decir, la tasa de disminución queda tal que:
$$ \frac{dy}{dt} = -\gamma y_{0} e^{-\gamma t} \;\;\;\Rightarrow\;\;\; \frac{dy}{dt} = -\gamma y $$
Al igual que con las presas, considerando un análisis conjunto. Tenemos que la población de depredadores se ve afectada por la población de presas, permitiendo su proliferación gracias a la disponibilidad de alimento. Este efecto relaciona las poblaciones de presas y depredadores mediante una eficiencia de conversión de presas en nuevos depredadores $(\delta)$. Finalmente, la tasa de crecimiento de los depredadores se expresa como:
$$\frac{dy}{dt} = \delta xy - \gamma y$$
Para cada una de las constantes $\alpha$, $\beta$, $\gamma$ y $\delta$, existen formas prácticas de estimación. La primera consiste en realizar un experimento controlado en el que, en primer lugar, se analiza cada población de manera aislada para determinar $\alpha$ y $\gamma$, y posteriormente se estudia el comportamiento conjunto de ambas poblaciones para estimar $\beta$ y $\delta$. Otra manera de aproximar los valores de las constantes es obtener las series de tiempo $x(t)$ e $y(t)$ de una población real y ajustar el modelo a los datos, de modo de obtener los parámetros que mejor los expliquen. De esta forma, el modelo de Lotka-Volterra queda como $$ \begin{cases} \frac{dx}{dt} &= \alpha x - \beta xy \\ \frac{dy}{dt} &= \delta xy - \gamma y \end{cases} $$
A partir de este punto utilizaremos la siguiente notación para el modelo Lotka–Volterra:
- $p(t)$: población de presas en el tiempo $t$
- $P(t)$: población de depredadores en el tiempo $t$
Con sus respectivas constantes:
- $\alpha$: tasa de crecimiento natural de las presas.
- $\beta$: tasa de depredación (probabilidad de encuentro entre presas y depredadores).
- $\epsilon$: eficiencia de conversión de presas en nuevos depredadores.
- $m$: tasa de mortalidad natural de los depredadores.
El sistema de ecuaciones diferenciales queda expresado como: $$ \begin{cases} \frac{dp}{dt} &= \alpha p - \beta p P \\ \frac{dP}{dt} &= \epsilon \beta p P - m P \end{cases}
1. Implementación Numérica¶
# Función para modelar la tasa de población de las presas.
def prey_poblation(a: float, b: float, p: int, P: int) -> float:
"""
Calcula la tasa de cambio de la población de presas según el modelo Lotka-Volterra.
Parámetros:
----------
a (float): Tasa de crecimiento de las presas.
b (float): Tasa de depredación.
p (int): Población de presas.
P (int): Población de depredadores.
Retorna:
-------
float: Tasa de cambio de la población de presas.
"""
return (a * p) - (b * p * P)
# Función para modelar la tasa de población de los depredadores.
def depr_poblation(b: float, m: float, p: int, P: int, eps: float) -> float:
"""
Calcula la tasa de cambio de la población de depredadores según el modelo Lotka-Volterra.
Parámetros:
----------
b (float): Tasa de depredación.
m (float): Tasa de mortalidad de depredadores.
p (int): Población de presas.
P (int): Población de depredadores.
eps (float): Eficiencia de conversión de presas en depredadores.
Retorna:
-------
float: Tasa de cambio de la población de depredadores.
"""
return (eps * b * p * P) - (m * P)
# Función para modelar la dinámica del sistema Lotka-Volterra.
def lotka_volterra(t: float, y: list[float], params: dict) -> list[float]:
"""
Calcula las derivadas de las poblaciones de presas y depredadores en el modelo Lotka-Volterra.
Parámetros:
----------
t (float): Tiempo actual (no se usa en el cálculo).
y (list[float]): Lista con las poblaciones actuales [presas, depredadores].
params (dict): Diccionario con los parámetros del modelo (a, b, m, eps).
Retorna:
-------
list[float]: Lista con las tasas de cambio [dp/dt, dP/dt].
"""
p, P = y
a, b, m, eps = params['a'], params['b'], params['m'], params['eps']
dpdt = prey_poblation(a, b, p, P)
dPdt = depr_poblation(b, m, p, P, eps)
return [dpdt, dPdt]
def simulate(params: dict, y0: list[float], dt: float, T: int, method: Callable) -> dict[str, list]:
"""
Simula el sistema Lotka-Volterra utilizando un método numérico.
Parámetros:
----------
params (dict): Diccionario con los parámetros del modelo (a, b, m, eps).
y0 (list[float]): Condiciones iniciales [presas, depredadores].
dt (float): Paso de tiempo.
T (int): Tiempo total de simulación.
method (Callable): Función que implementa el método numérico (e.g., euler, rk2, rk4).
Retorna:
-------
dict[str, list]: Diccionario con los resultados de la simulación:
'times': lista de tiempos,
'preys': lista de poblaciones de presas,
'predators': lista de poblaciones de depredadores.
"""
steps = int(np.floor(T/dt))
# Almacenamos las condiciones iniciales.
y = y0
# Listas para almacenar los resultados.
times = []
preys = []
predators = []
t = 0
for _ in range(steps):
# Almacenamos los tiempos.
times.append(t)
# Almacenamos las poblaciones de presas para cada método.
preys.append(y[0])
# Almacenamos las poblaciones de depredadores para cada método.
predators.append(y[1])
# Actualizamos las poblaciones usando cada método numérico.
y = method(lotka_volterra, t, y, dt, params)
t += dt
# Almacenamos los resultados
results = {
'times': times,
'preys': preys,
'predators': predators
}
return results
Funciones para los métodos numéricos¶
# Implementación del método de Euler.
def euler_step(f: Callable, t: float, y: list[float], dt: float, params: dict) -> list[float]:
"""
Realiza un paso del método de Euler para resolver una ecuación diferencial ordinaria.
Parámetros:
----------
f (Callable): Función que calcula las derivadas.
t (float): Tiempo actual.
y (list[float]): Lista con los valores actuales de las variables dependientes.
dt (float): Paso de tiempo.
params (dict): Diccionario con parámetros adicionales para la función f.
Retorna:
-------
list[float]: Lista con los nuevos valores de las variables dependientes después del paso.
"""
dydt = f(t, y, params)
return [y[i] + dt * dydt[i] for i in range(len(y))]
# Implementación del método de Runge-Kutta de segundo orden (RK2).
def rk2_step(f: Callable, t: float, y: list[float], dt: float, params: dict) -> list[float]:
"""
Realiza un paso del método de Runge-Kutta de segundo orden para resolver una ecuación diferencial ordinaria.
Parámetros:
----------
f (Callable): Función que calcula las derivadas.
t (float): Tiempo actual.
y (list[float]): Lista con los valores actuales de las variables dependientes.
dt (float): Paso de tiempo.
params (dict): Diccionario con parámetros adicionales para la función f.
Retorna:
-------
list[float]: Lista con los nuevos valores de las variables dependientes después del paso.
"""
k1 = [dt * val for val in f(t, y, params)]
k2 = [dt * val for val in f(t + dt, [y[i] + k1[i] for i in range(len(y))], params)]
return [y[i] + (k1[i] + k2[i]) / 2 for i in range(len(y))]
# Implementación del método de Runge-Kutta de cuarto orden (RK4).
def rk4_step(f: Callable, t: float, y: list[float], dt: float, params: dict) -> list[float]:
"""
Realiza un paso del método de Runge-Kutta de cuarto orden para resolver una ecuación diferencial ordinaria.
Parámetros:
----------
f (Callable): Función que calcula las derivadas.
t (float): Tiempo actual.
y (list[float]): Lista con los valores actuales de las variables dependientes.
dt (float): Paso de tiempo.
params (dict): Diccionario con parámetros adicionales para la función f.
Retorna:
-------
list[float]: Lista con los nuevos valores de las variables dependientes después del paso.
"""
k1 = [dt * val for val in f(t, y, params)]
k2 = [dt * val for val in f(t + (dt/2), [y[i] + (k1[i]/2) for i in range(len(y))], params)]
k3 = [dt * val for val in f(t + (dt/2), [y[i] + (k2[i]/2) for i in range(len(y))], params)]
k4 = [dt * val for val in f(t + dt, [y[i] + k3[i]for i in range(len(y))], params)]
return [y[i] + ((k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) / 6) for i in range(len(y))]
Comparación de las soluciones de los diferentes métodos numéricos.¶
Simulación del modelo Lotka-Volterra con diferentes métodos numéricos¶
# Parámetros del modelo Lotka-Volterra. Fueron elegidos de forma arbitraria.
params = {
'a': 0.05,
'b': 0.003,
'm': 0.01,
'eps': 0.02,
}
# Condiciones iniciales: 40 presas y 25 depredadores.
y0 = [40, 25]
# Configuración de la simulación. Paso de tiempo y número total de pasos.
dt = 0.01
T = 50_000
results_eu = simulate(params, y0, dt, T, euler_step)
results_rk2 = simulate(params, y0, dt, T, rk2_step)
results_rk4 = simulate(params, y0, dt, T, rk4_step)
# Cálculo de las relaciones de equilibrio.
R_e = params['m'] / (params['eps'] * params['b'])
F_e = params['a'] / params['b']
Visualización de los resultados¶
plt.figure(figsize=(16, 18))
plt.subplot(3, 1, 1)
plt.plot(results_eu['times'], results_eu['preys'], label='Presas (Euler)', color='blue', alpha=0.5)
plt.plot(results_eu['times'], results_eu['predators'], label='Depredadores (Euler)', color='red', alpha=0.5)
plt.xlabel('Tiempo')
plt.xlim(0, 2_000)
plt.ylabel('Población')
plt.legend(loc='upper right')
plt.title('Método de Euler')
plt.subplot(3, 1, 2)
plt.plot(results_rk2['times'], results_rk2['preys'], label='Presas (RK2)', color='blue', alpha=0.5)
plt.plot(results_rk2['times'], results_rk2['predators'], label='Depredadores (RK2)', color='red', alpha=0.5)
plt.xlabel('Tiempo')
plt.xlim(0, 2_000)
plt.ylabel('Población')
plt.legend(loc='upper right')
plt.title('Método de Runge-Kutta 2')
plt.subplot(3, 1, 3)
plt.plot(results_rk4['times'], results_rk4['preys'], label='Presas (RK4)', color='blue', alpha=0.5)
plt.plot(results_rk4['times'], results_rk4['predators'], label='Depredadores (RK4)', color='red', alpha=0.5)
plt.xlabel('Tiempo')
plt.xlim(0, 2_000)
plt.ylabel('Población')
plt.legend(loc='upper right')
plt.title('Método de Runge-Kutta 4')
plt.tight_layout()
plt.show();
plt.figure(figsize=(16, 16))
plt.subplot(2, 2, 1)
plt.plot(results_eu['times'], results_eu['preys'], label='Presas (Euler)', color='blue', alpha=0.5)
plt.plot(results_rk2['times'], results_rk2['preys'], label='Presas (RK2)', color='green', alpha=0.5)
plt.plot(results_rk4['times'], results_rk4['preys'], label='Presas (RK4)', color='magenta', alpha=0.5)
plt.xlabel('Tiempo')
plt.xlim(0, 2_000)
plt.ylabel('Población de presas')
plt.legend(loc='upper right')
plt.title('Comparación de métodos numéricos en un delta de tiempo grande\npara la población de presas.')
plt.subplot(2, 2, 2)
plt.plot(results_eu['times'], results_eu['predators'], label='Depredadores (Euler)', color='blue', alpha=0.5)
plt.plot(results_rk2['times'], results_rk2['predators'], label='Depredadores (RK2)', color='green', alpha=0.5)
plt.plot(results_rk4['times'], results_rk4['predators'], label='Depredadores (RK4)', color='magenta', alpha=0.5)
plt.xlabel('Tiempo')
plt.xlim(0, 2_000)
plt.ylabel('Población de depredadores')
plt.legend(loc='upper right')
plt.title('Comparación de métodos numéricos en un delta de tiempo grande\npara la población de depredadores.')
plt.subplot(2, 2, 3)
plt.plot(results_eu['times'], results_eu['preys'], label='Presas (Euler)', color='orange', alpha=0.7)
plt.plot(results_rk2['times'], results_rk2['preys'], label='Presas (RK2)', color='red', alpha=0.7)
plt.plot(results_rk4['times'], results_rk4['preys'], label='Presas (RK4)', color='darkorange', alpha=0.7)
plt.xlabel('Tiempo')
plt.xlim(1_790, 1_800)
plt.ylabel('Población de presas')
plt.ylim(540, 565)
plt.yticks(np.arange(540, 566, 5))
plt.legend(loc='upper right')
plt.title('Comparación de métodos numéricos en un delta de tiempo pequeño\npara la población de presas.')
plt.subplot(2, 2, 4)
plt.plot(results_eu['times'], results_eu['predators'], label='Depredadores (Euler)', color='cyan', alpha=0.7)
plt.plot(results_rk2['times'], results_rk2['predators'], label='Depredadores (RK2)', color='deepskyblue', alpha=0.7)
plt.plot(results_rk4['times'], results_rk4['predators'], label='Depredadores (RK4)', color='navy', alpha=0.7)
plt.xlabel('Tiempo')
plt.xlim(1_790, 1_800)
plt.ylabel('Población de depredadores')
plt.ylim(10, 20)
plt.legend(loc='upper right')
plt.title('Comparación de métodos numéricos en un delta de tiempo pequeño\npara la población de depredadores.')
plt.tight_layout()
plt.show();
En las sumulaciones se compararon 3 métodos de integración numérica: Euler, Runge-Kutta de orden 2 $(RK2)$ y Runge-Kutta de orden 4 $(RK4)$ en un intervalo de tiempo grande y uno pequeño, con el fin de obtener una visión macro y micro de las curvas en función del tiempo.
Para un intervalo de tiempo grande se observa que los 3 métodos producen oscilaciones periódicas en las poblaciones de presas y depredadores, en donde todos los métodos describen el mismo comportamiento, con variaciones que a esa escala no se perciben. Las curvas, a simple vista, se superponen mostrando trayectorias similares, esto es a pesar de la diferencia en la presición de los métodos, dado que incluso el método con menos presición (Euler) logra capturar satisfactoriamente la dinámica general del sistema.
En un intervalo de tiempo más pequeño, para una mirada más local de las curvas, las variaciones y diferencias entre los métodos se hacen más evidentes. A una escala local, el método de Euler suele sobreestimar o subestimar en comparación con los métodos de Runge–Kutta. Respecto a estos dos últimos, a pesar de ser de distinto orden, presentan una superposición que los hace prácticamente indistinguibles en ese rango de valores de $t$.
Funciones auxiliares para la animación.¶
# Con la ayuda de ChatGPT con pensamiento.
def simulate_for_animation(params: dict, p0: float, P0: float, dt: float, T_anim: float, frames_target: int = 1000, clip_negative: bool = False) -> dict[str, np.ndarray]:
steps = int(np.ceil(T_anim / dt))
save_every = max(1, steps // frames_target)
times = []
prey_eu = []; prey_rk2 = []; prey_rk4 = []
pred_eu = []; pred_rk2 = []; pred_rk4 = []
y_eu = [p0, P0]; y_rk2 = [p0, P0]; y_rk4 = [p0, P0]
t = 0.0
for step in range(steps + 1):
if step % save_every == 0:
times.append(t)
prey_eu.append(y_eu[0]); prey_rk2.append(y_rk2[0]); prey_rk4.append(y_rk4[0])
pred_eu.append(y_eu[1]); pred_rk2.append(y_rk2[1]); pred_rk4.append(y_rk4[1])
# un paso
y_eu = euler_step(lotka_volterra, t, y_eu, dt, params)
y_rk2 = rk2_step(lotka_volterra, t, y_rk2, dt, params)
y_rk4 = rk4_step(lotka_volterra, t, y_rk4, dt, params)
if clip_negative:
y_eu = [max(0.0, v) for v in y_eu]
y_rk2 = [max(0.0, v) for v in y_rk2]
y_rk4 = [max(0.0, v) for v in y_rk4]
t += dt
return {
'times': np.array(times),
'prey_eu': np.array(prey_eu),'prey_rk2': np.array(prey_rk2),'prey_rk4': np.array(prey_rk4),
'pred_eu': np.array(pred_eu),'pred_rk2': np.array(pred_rk2),'pred_rk4': np.array(pred_rk4),
}
def make_animation_prey(params, p0, P0, dt=0.01, T_anim=200.0, frames_target=600, clip_negative=False, fps=15, zoom_window=20.0, figsize=(10,4)) -> FuncAnimation:
data = simulate_for_animation(params, p0, P0, dt, T_anim, frames_target, clip_negative)
times = data['times']
N = len(times)
fig, (ax_full, ax_zoom) = plt.subplots(1,2, figsize=figsize)
fig.suptitle(f'Presas: Euler / RK2 / RK4 (dt={dt}, T={T_anim})')
# --- full plot (serie completa) ---
line_full_eu, = ax_full.plot([], [], lw=1, alpha=0.6, label='Euler')
line_full_rk2, = ax_full.plot([], [], lw=1, alpha=0.8, label='RK2')
line_full_rk4, = ax_full.plot([], [], lw=1, alpha=0.9, label='RK4')
mark_full_eu, = ax_full.plot([], [], marker='o', ms=5, alpha=0.6)
mark_full_rk2, = ax_full.plot([], [], marker='o', ms=5, alpha=0.8)
mark_full_rk4, = ax_full.plot([], [], marker='o', ms=5, alpha=0.9)
ax_full.set_xlabel('Tiempo'); ax_full.set_ylabel('Población de presas')
ax_full.legend(loc='upper right')
ax_full.set_xlim(times[0], times[-1])
# y-limits según datos (añadimos pequeño padding)
ymin = min(data['prey_eu'].min(), data['prey_rk2'].min(), data['prey_rk4'].min())
ymax = max(data['prey_eu'].max(), data['prey_rk2'].max(), data['prey_rk4'].max())
ax_full.set_ylim(max(0, ymin*0.9), ymax*1.1)
# --- zoom plot (ventana móvil en tiempo) ---
line_zoom_eu, = ax_zoom.plot([], [], lw=1, alpha=0.6, label='Euler')
line_zoom_rk2, = ax_zoom.plot([], [], lw=1, alpha=0.8, label='RK2')
line_zoom_rk4, = ax_zoom.plot([], [], lw=1, alpha=0.9, label='RK4')
mark_zoom_eu, = ax_zoom.plot([], [], marker='o', ms=6, alpha=0.6)
mark_zoom_rk2, = ax_zoom.plot([], [], marker='o', ms=6, alpha=0.8)
mark_zoom_rk4, = ax_zoom.plot([], [], marker='o', ms=6, alpha=0.9)
ax_zoom.set_xlabel('Tiempo'); ax_zoom.set_ylabel('Población (zoom)')
ax_zoom.legend(loc='upper right')
def init():
for ln in [line_full_eu, line_full_rk2, line_full_rk4,
mark_full_eu, mark_full_rk2, mark_full_rk4,
line_zoom_eu, line_zoom_rk2, line_zoom_rk4,
mark_zoom_eu, mark_zoom_rk2, mark_zoom_rk4]:
ln.set_data([], [])
return (line_full_eu, line_full_rk2, line_full_rk4,
mark_full_eu, mark_full_rk2, mark_full_rk4,
line_zoom_eu, line_zoom_rk2, line_zoom_rk4,
mark_zoom_eu, mark_zoom_rk2, mark_zoom_rk4)
def update(frame):
# frame corresponde al índice en arrays muestreados
t_all = times
t_current = times[frame]
# series hasta el frame actual (para la gráfica completa)
t = t_all[:frame+1]
p_eu = data['prey_eu'][:frame+1]
p_rk2 = data['prey_rk2'][:frame+1]
p_rk4 = data['prey_rk4'][:frame+1]
# actualizar full (izquierda)
line_full_eu.set_data(t, p_eu)
line_full_rk2.set_data(t, p_rk2)
line_full_rk4.set_data(t, p_rk4)
mark_full_eu.set_data([t_current], [data['prey_eu'][frame]])
mark_full_rk2.set_data([t_current], [data['prey_rk2'][frame]])
mark_full_rk4.set_data([t_current], [data['prey_rk4'][frame]])
# --- Ventana móvil centrada en t_current (derecha) ---
# Define ventana centrada
half_win = zoom_window / 2.0
t0 = t_current - half_win
t1 = t_current + half_win
# clamp a límites de los datos
t_min_all, t_max_all = t_all[0], t_all[-1]
if t0 < t_min_all:
t0 = t_min_all
t1 = t_min_all + zoom_window
if t1 > t_max_all:
t1 = t_max_all
t0 = max(t_min_all, t_max_all - zoom_window)
# mascara para datos dentro de la ventana centrada
t_mask = (t_all >= t0) & (t_all <= t1)
times_zoom = t_all[t_mask]
# datos correspondientes al zoom (podrían estar más espaciados si muestreaste)
p_eu_zoom = data['prey_eu'][t_mask]
p_rk2_zoom = data['prey_rk2'][t_mask]
p_rk4_zoom = data['prey_rk4'][t_mask]
# actualizar líneas del zoom
line_zoom_eu.set_data(times_zoom, p_eu_zoom)
line_zoom_rk2.set_data(times_zoom, p_rk2_zoom)
line_zoom_rk4.set_data(times_zoom, p_rk4_zoom)
mark_zoom_eu.set_data([t_current], [data['prey_eu'][frame]])
mark_zoom_rk2.set_data([t_current], [data['prey_rk2'][frame]])
mark_zoom_rk4.set_data([t_current], [data['prey_rk4'][frame]])
# --- Ajuste dinámico de límites Y para centrar el punto ---
# Parámetros que puedes ajustar (si los mueves a la firma, podrás controlarlos desde fuera)
zoom_y_factor = 0.6 # 0.6 -> la mitad (60%) del rango local; menor => zoom más apretado verticalmente
min_zoom_y = 0.5 # límite mínimo para el semi-rango Y (evita 0)
if len(p_eu_zoom) > 0:
# rango local observado
local_y_min = min(p_eu_zoom.min(), p_rk2_zoom.min(), p_rk4_zoom.min())
local_y_max = max(p_eu_zoom.max(), p_rk2_zoom.max(), p_rk4_zoom.max())
local_range = max(1e-8, local_y_max - local_y_min)
# centro vertical: puedes usar la media local o el valor actual. Aquí usamos el valor actual promedio entre métodos
cur_y_vals = np.array([data['prey_eu'][frame], data['prey_rk2'][frame], data['prey_rk4'][frame]])
y_center = float(np.mean(cur_y_vals))
# semi-rango vertical deseado: proporcional al rango local pero centrado en y_center
half_y = max(min_zoom_y, (local_range * zoom_y_factor) / 2.0)
# límites Y centrados en y_center (evitan negativos)
y0 = max(0.0, y_center - half_y)
y1 = y_center + half_y
# si la ventana vertical queda demasiado pequeña en relación a los datos (p. ej. todos iguales), ampliar un poco
if (y1 - y0) < 1e-3:
y0 = max(0.0, y_center - min_zoom_y)
y1 = y_center + min_zoom_y
ax_zoom.set_xlim(t0, t1)
ax_zoom.set_ylim(y0, y1)
else:
# fallback si no hay datos en la ventana (raro)
ax_zoom.set_xlim(t0, t1)
return (line_full_eu, line_full_rk2, line_full_rk4,
mark_full_eu, mark_full_rk2, mark_full_rk4,
line_zoom_eu, line_zoom_rk2, line_zoom_rk4,
mark_zoom_eu, mark_zoom_rk2, mark_zoom_rk4)
ani = animation.FuncAnimation(fig, update, frames=N, init_func=init, blit=False, interval=1000/fps)
plt.tight_layout()
return ani
Visualización de la animación.¶
params = {
'a': 0.05,
'b': 0.003,
'm': 0.01,
'eps': 0.02,
}
p0, P0 = 40.0, 25.0
# Parámetros para la animación
dt = 0.1
T_anim = 400.0
frames_target = 480
fps = 24
bitrate = 4_000
dpi = 150
zoom_window = 30.0
clip_negative = True
# Generar la animación
ani = make_animation_prey(params, p0, P0, dt=dt, T_anim=T_anim,
frames_target=frames_target, clip_negative=clip_negative,
fps=fps, zoom_window=zoom_window, figsize=(12, 6))
# Guardar la animación en un archivo MP4
part1_animation_fname = 'media/prey_zoom.mp4' # Nombre del archivo de salida
# Verifica que ffmpeg esté disponible
ff = shutil.which('ffmpeg')
if ff is None:
raise RuntimeError("ffmpeg no encontrado en PATH. " \
"Ejecuta 'where ffmpeg' en CMD o añade la ruta a mpl.rcParams['animation.ffmpeg_path'].")
# Guardar usando FFMpegWriter
writer = animation.FFMpegWriter(fps=fps, bitrate=bitrate)
ani.save(part1_animation_fname, writer=writer, dpi=dpi)
size = os.path.getsize(part1_animation_fname) / (1024*1024)
print(f"Animación guardada en '{part1_animation_fname}' con tamaño {size:.2f} MB.")
plt.close()
Animación guardada en 'media/prey_zoom.mp4' con tamaño 5.92 MB.
display(Video('media/prey_zoom.mp4', embed=True, width=1200))
2. Análisis del Paso de Tiempo¶
def simulate_compare(params: dict, p0: float, P0: float, dt_list: list[float], T: float, save_max_points: int = 10000, clip_negative: bool = False):
results = {}
for dt_i in dt_list:
steps = int(T / dt_i)
# advertencia si demasiados pasos
if steps > 5_000_000:
print(f"Warning: dt={dt_i} -> steps={steps} (muy grande). Considera reducir T o aumentar dt.")
# decidir cada cuantos pasos guardo (muestreo)
save_every = max(1, steps // save_max_points)
# reiniciar condiciones para este dt
y_eu = [p0, P0]
y_rk2 = [p0, P0]
y_rk4 = [p0, P0]
times = []
prey_eu, prey_rk2, prey_rk4 = [], [], []
pred_eu, pred_rk2, pred_rk4 = [], [], []
t = 0.0
for step in range(steps):
# guardar cada save_every pasos (incluyendo t=0)
if step % save_every == 0:
times.append(t)
prey_eu.append(y_eu[0])
prey_rk2.append(y_rk2[0])
prey_rk4.append(y_rk4[0])
pred_eu.append(y_eu[1])
pred_rk2.append(y_rk2[1])
pred_rk4.append(y_rk4[1])
# dar un paso
y_eu = euler_step(lotka_volterra, t, y_eu, dt_i, params)
y_rk2 = rk2_step(lotka_volterra, t, y_rk2, dt_i, params)
y_rk4 = rk4_step(lotka_volterra, t, y_rk4, dt_i, params)
if clip_negative:
y_eu = [max(0.0, v) for v in y_eu]
y_rk2 = [max(0.0, v) for v in y_rk2]
y_rk4 = [max(0.0, v) for v in y_rk4]
t += dt_i
results[f"dt_{dt_i}"] = {
'times': times,
'prey_eu': prey_eu, 'prey_rk2': prey_rk2, 'prey_rk4': prey_rk4,
'pred_eu': pred_eu, 'pred_rk2': pred_rk2, 'pred_rk4': pred_rk4,
'save_every': save_every,
'steps': steps
}
return results
params = {
'a': 0.05,
'b': 0.003,
'm': 0.01,
'eps': 0.02,
}
p0, P0 = 40, 25
dt = [0.5, 0.1, 0.01, 0.001]
T = 5_000
results = simulate_compare(params, p0, P0, dt, T, save_max_points=2000, clip_negative=True)
n = len(dt)
fig, axes = plt.subplots(nrows=n, ncols=2, figsize=(16, 5*n), sharex=False)
fig.suptitle(r'Comparación: Euler vs RK2 vs RK4 para distintos $\Delta t$', fontsize=16)
for i, dt_i in enumerate(dt):
res = results[f"dt_{dt_i}"]
times = res['times']
prey_eu = res['prey_eu']
prey_rk2 = res['prey_rk2']
prey_rk4 = res['prey_rk4']
pred_eu = res['pred_eu']
pred_rk2 = res['pred_rk2']
pred_rk4 = res['pred_rk4']
ax_prey = axes[i, 0] if n > 1 else axes[0]
ax_pred = axes[i, 1] if n > 1 else axes[1]
ax_prey.plot(times, prey_eu, label='Presas (Euler)', alpha=0.6)
ax_prey.plot(times, prey_rk2, label='Presas (RK2)', alpha=0.7)
ax_prey.plot(times, prey_rk4, label='Presas (RK4)', alpha=0.9)
ax_prey.set_ylabel('Población de presas')
ax_prey.set_title(f'Presas (dt={dt_i})')
ax_prey.legend(loc='upper left')
ax_pred.plot(times, pred_eu, label='Depredadores (Euler)', alpha=0.6)
ax_pred.plot(times, pred_rk2, label='Depredadores (RK2)', alpha=0.7)
ax_pred.plot(times, pred_rk4, label='Depredadores (RK4)', alpha=0.9)
ax_pred.set_ylabel('Población de depredadores')
ax_pred.set_title(f'Depredadores (dt={dt_i})')
ax_pred.legend(loc='upper left')
for ax in fig.axes:
ax.set_xlabel('Tiempo')
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()
En el análisis de la influencia de $\Delta t$, se observa que su elección impacta directamente en la estabilidad y fidelidad de las soluciones obtenidas. Se pueden observar distintas curvas asociadas a ambas poblaciones para los valores 0.5, 0.1, 0.01 y 0.001 de $\Delta t$.
Para $\Delta t = 0.5$, el método de Euler comienza a desviarse notoriamente respecto a los métodos de Runge-Kutta. Las desviaciones generan oscilaciones con amplitudes no realistas y acumulación de error a lo largo del tiempo. Por otro lado, $RK2$ y $RK4$ mantienen trayectorias mucho más estables y cercanas entre sí, capturando la dinámica poblacional esperada del sistema.
Para $\Delta t = 0.1$ las diferencias entre lo métodos disminuyen, aunque aún se aprecian discrepancias entre Euler y los métodos de mayor orden
Para $\Delta t$ más pequeños, como 0.01 y 0.001, las soluciones obtenidas con los 3 métodos son casi indistinguibles y describen practicamente las mismas curvas de poblaciones. A esta escala el método de Euler logra describir de manera adecuada la dinámica esperada del modelo, aunque con menos precisión que lo modelos de Runge-Kutta, tal como se mencionó en el inciso anterior (En la que se utilizó un $\Delta t = 0.01$).
3. Análisis de Error¶
def simulate_single(method_step: Callable, params: dict, y0: List[float],
dt: float, T: float, clip_negative: bool=False) -> tuple[np.ndarray, np.ndarray]:
steps = int(np.floor(T / dt))
times = np.linspace(0.0, steps * dt, steps + 1)
y = list(y0)
traj = np.zeros((steps + 1, 2))
traj[0, :] = y
t = 0.0
for i in range(1, steps+1):
y = method_step(lotka_volterra, t, y, dt, params)
if clip_negative:
y = [max(0.0, v) for v in y]
t += dt
traj[i, :] = y
return times, traj
def instant_error(y_ref: np.ndarray, y_method: np.ndarray) -> np.ndarray:
diff = y_ref - y_method
l2 = np.sqrt(np.sum(diff**2, axis=1))
return l2
def compare_methods(params: dict, y0: list[float], dt: float, T: float, clip_negative: bool=False,
solver_method: str='DOP853', rtol: float=1e-13, atol: float=1e-15) -> dict[str, any]:
times, traj_eu = simulate_single(euler_step, params, y0, dt, T, clip_negative=clip_negative)
_, traj_rk2 = simulate_single(rk2_step, params, y0, dt, T, clip_negative=clip_negative)
_, traj_rk4 = simulate_single(rk4_step, params, y0, dt, T, clip_negative=clip_negative)
t0 = 0.0
tf = times[-1]
def f_wrap(t: float, y: np.ndarray) -> np.ndarray:
return np.array(lotka_volterra(t, list(y), params), dtype=float)
sol = solve_ivp(fun=f_wrap, t_span=(0, T), y0=y0, method=solver_method,
t_eval=times, rtol=rtol, atol=atol)
if not sol.success:
raise RuntimeError("solve_ivp falló: " + str(sol.message))
traj_ref = sol.y.T
err_eu = instant_error(traj_ref, traj_eu)
err_rk2 = instant_error(traj_ref, traj_rk2)
err_rk4 = instant_error(traj_ref, traj_rk4)
stats = {
'dt': dt,
'times': times,
'traj_ref': traj_ref,
'traj_eu': traj_eu, 'traj_rk2': traj_rk2, 'traj_rk4': traj_rk4,
'err_eu': err_eu, 'err_rk2': err_rk2, 'err_rk4': err_rk4,
'mean_err': (np.mean(err_eu), np.mean(err_rk2), np.mean(err_rk4)),
'max_err': (np.max(err_eu), np.max(err_rk2), np.max(err_rk4)),
'rms_err': (np.sqrt(np.mean(err_eu**2)), np.sqrt(np.mean(err_rk2**2)), np.sqrt(np.mean(err_rk4**2))),
}
return stats
params = {
'a': 0.05,
'b': 0.003,
'm': 0.01,
'eps': 0.02,
'K': 1_500
}
p0, P0 = 40, 25
y0 = [p0, P0]
T = 200
dt_list = [0.5, 0.1, 0.05, 0.01, 0.005, 0.001, 0.0005, 0.0001]
results_err = []
for dt in dt_list:
stats = compare_methods(params, y0, dt, T, clip_negative=False)
results_err.append(stats)
print('dt:', dt)
print('Error cuadrático medio (Euler, RK2, RK4):', np.round(stats['rms_err'], 16))
if dt != dt_list[-1]:
print()
dt: 0.5 Error cuadrático medio (Euler, RK2, RK4): [1.27711468e+00 7.76953133e-04 3.72545224e-08] dt: 0.1 Error cuadrático medio (Euler, RK2, RK4): [2.52009017e-01 3.15233002e-05 7.67524000e-11] dt: 0.05 Error cuadrático medio (Euler, RK2, RK4): [1.25790889e-01 7.89518733e-06 3.80854000e-11] dt: 0.01 Error cuadrático medio (Euler, RK2, RK4): [2.51239920e-02 3.16270402e-07 3.64984000e-11] dt: 0.005 Error cuadrático medio (Euler, RK2, RK4): [1.25598592e-02 7.90823008e-08 3.77637000e-11] dt: 0.001 Error cuadrático medio (Euler, RK2, RK4): [2.51162995e-03 3.16621140e-09 3.67580000e-11] dt: 0.0005 Error cuadrático medio (Euler, RK2, RK4): [1.2557936e-03 7.8485460e-10 3.7047500e-11] dt: 0.0001 Error cuadrático medio (Euler, RK2, RK4): [2.51155327e-04 5.94650000e-11 4.07932000e-11]
dt_vals = [res['dt'] for res in results_err]
rms_eu = [res['rms_err'][0] for res in results_err]
rms_rk2 = [res['rms_err'][1] for res in results_err]
rms_rk4 = [res['rms_err'][2] for res in results_err]
plt.figure(figsize=(12, 6))
plt.plot(dt_vals, rms_eu, 'o-', label='Euler (RMS)')
plt.plot(dt_vals, rms_rk2, 's-', label='RK2 (RMS)')
plt.plot(dt_vals, rms_rk4, '^-', label='RK4 (RMS)')
plt.gca().invert_xaxis()
plt.xlabel('Paso de Tiempo (dt)')
plt.xscale('log')
plt.ylabel('Error RMS')
plt.yscale('log')
plt.title('Convergencia de Métodos Numéricos')
plt.legend()
plt.grid(True, which="both", linewidth=0.5)
plt.show();
def local_orders(dt_vals, errs):
orders = []
for i in range(len(dt_vals)-1):
h1, h2 = dt_vals[i], dt_vals[i+1]
e1, e2 = errs[i], errs[i+1]
if e1>0 and e2>0:
p = np.log(e1/e2)/np.log(h1/h2)
else:
p = np.nan
orders.append(p)
return np.array(orders)
orders_eu = local_orders(dt_vals, rms_eu)
orders_rk2 = local_orders(dt_vals, rms_rk2)
orders_rk4 = local_orders(dt_vals, rms_rk4)
for name, ords in [('Euler', orders_eu), ('RK2', orders_rk2), ('RK4', orders_rk4)]:
print(f"Órdenes locales de convergencia para {name}:")
for i, p in enumerate(ords):
print(f"-> Entre dt={dt_vals[i]} y dt={dt_vals[i+1]}: p ≈ {p:.2f}")
if name != 'RK4':
print()
Órdenes locales de convergencia para Euler: -> Entre dt=0.5 y dt=0.1: p ≈ 1.01 -> Entre dt=0.1 y dt=0.05: p ≈ 1.00 -> Entre dt=0.05 y dt=0.01: p ≈ 1.00 -> Entre dt=0.01 y dt=0.005: p ≈ 1.00 -> Entre dt=0.005 y dt=0.001: p ≈ 1.00 -> Entre dt=0.001 y dt=0.0005: p ≈ 1.00 -> Entre dt=0.0005 y dt=0.0001: p ≈ 1.00 Órdenes locales de convergencia para RK2: -> Entre dt=0.5 y dt=0.1: p ≈ 1.99 -> Entre dt=0.1 y dt=0.05: p ≈ 2.00 -> Entre dt=0.05 y dt=0.01: p ≈ 2.00 -> Entre dt=0.01 y dt=0.005: p ≈ 2.00 -> Entre dt=0.005 y dt=0.001: p ≈ 2.00 -> Entre dt=0.001 y dt=0.0005: p ≈ 2.01 -> Entre dt=0.0005 y dt=0.0001: p ≈ 1.60 Órdenes locales de convergencia para RK4: -> Entre dt=0.5 y dt=0.1: p ≈ 3.84 -> Entre dt=0.1 y dt=0.05: p ≈ 1.01 -> Entre dt=0.05 y dt=0.01: p ≈ 0.03 -> Entre dt=0.01 y dt=0.005: p ≈ -0.05 -> Entre dt=0.005 y dt=0.001: p ≈ 0.02 -> Entre dt=0.001 y dt=0.0005: p ≈ -0.01 -> Entre dt=0.0005 y dt=0.0001: p ≈ -0.06
El análisis de error se realizó comparando las soluciones obtenidas mediante los métodos de Euler, $RK2$ y $RK4$ respecto a una solución de referencia generada por un integrador de mayor orden, específicamente, un método Runge-Kutta 8 (DOP853) con tolerancias estrictas, debido a que el modelo Lotka-Volterra no posee solución analítica, por tanto utilizamos el método anteriormente mencionado que proviene de la función solve_ivp de la librería scipy. Luego, la solución de referencia se evaluó en los mismos instantes de tiempo que las soluciones numéricas y se calculó el error mediante la norma RMS (raíz del error cuadrático medio). Por último se realizó un gráfico log-log del error en función de $\Delta t$ y las pendientes en diferentes puntos de cada curva, notando los errores característicos asociados a los modelos.
Para el método de Euler, se observa una pendiente $p \approx 1$ en el gráfico log–log de error, lo cual es coherente con su orden de convergencia lineal. En este caso, el error decrece proporcionalmente a $\Delta t$, por lo que al reducir el paso a la mitad, el error también se reduce aproximadamente a la mitad. Esto convierte a Euler en el método de menor precisión, ya que requiere pasos extremadamente pequeños para alcanzar errores aceptables.
Para el método Runge–Kutta de orden 2 (RK2), se observa una pendiente $p \approx 2$ en el gráfico log–log de error, lo que es consistente con su orden de convergencia cuadrático. En este caso, el error decrece proporcionalmente a $(\Delta t)^2$, de modo que si el paso se reduce a la mitad, el error se reduce aproximadamente a una cuarta parte. Esto lo hace considerablemente más preciso que el método de Euler, ya que permite obtener soluciones mucho más exactas sin necesidad de pasos tan pequeños, logrando un mejor equilibrio entre precisión y costo computacional.
Por último, para el método Runge–Kutta de orden 4 (RK4) se observa en primera instancia una pendiente $p \approx 4$ en el gráfico log–log de error, lo que corresponde a un decrecimiento proporcional a $(\Delta t)^4$. Esto significa que, si el paso de integración se reduce a la mitad, el error disminuye aproximadamente dieciséis veces. Este comportamiento confirma que RK4 es el método más preciso de los tres analizados, mostrando errores notablemente menores incluso con pasos de integración relativamente grandes y alcanzando una tasa de disminución mucho mayor que la del método de Euler y $RK2$ con $\Delta t$ más pequeños. Finalmente, la reducción del error se ve limitada por nuestra implementación del modelo poblacional, aunque esperariamos que si el modelo tuviese solución analítica, podríamos ver errores de precisión de máquinas con $\Delta t$ muy pequeños, esto explica la saturación observada en la parte baja de la curva ($\Delta t > 0.1$).
4. Diagrama de Fase¶
plt.figure(figsize=(16, 18))
plt.subplot(3, 2, 1)
plt.plot(results_eu['times'], results_eu['preys'], label='Presas (Euler)', color='blue')
plt.plot(results_eu['times'], results_eu['predators'], label='Depredadores (Euler)', color='red')
plt.xlabel('Tiempo')
plt.xlim(0, 2_000)
plt.ylabel('Población')
plt.legend(loc='upper right')
plt.title('Método de Euler')
plt.subplot(3, 2, 2)
plt.plot(results_eu['preys'], results_eu['predators'], color='purple')
plt.axvline(R_e, color='blue', linestyle='--', alpha=0.7, label='Equilibrio Presas')
plt.axhline(y=F_e, color='red', linestyle='--', alpha=0.7, label='Equilibrio Depredadores')
plt.scatter(R_e, F_e, color='black', zorder=5)
plt.annotate(f"Equilibrio:\nPresas={R_e:.1f}\nDepredadores={F_e:.1f}",
xy=(R_e, F_e), xytext=(R_e+80, F_e+7), fontsize=12,
arrowprops=dict(arrowstyle="->", color="black"))
plt.xlabel('Presas')
plt.ylabel('Depredadores')
plt.title('Espacio de fase con método de euler (presa v/s depredador)')
plt.subplot(3, 2, 3)
plt.plot(results_rk2['times'], results_rk2['preys'], label='Presas (RK2)', color='blue')
plt.plot(results_rk2['times'], results_rk2['predators'], label='Depredadores (RK2)', color='red')
plt.xlabel('Tiempo')
plt.xlim(0, 2_000)
plt.ylabel('Población')
plt.legend(loc='upper right')
plt.title('Método de Runge-Kutta 2')
plt.subplot(3, 2, 4)
plt.plot(results_rk2['preys'], results_rk2['predators'], color='purple')
plt.axvline(R_e, color='blue', linestyle='--', alpha=0.7, label='Equilibrio Presas')
plt.axhline(y=F_e, color='red', linestyle='--', alpha=0.7, label='Equilibrio Depredadores')
plt.scatter(R_e, F_e, color='black', zorder=5)
plt.annotate(f"Equilibrio:\nPresas={R_e:.1f}\nDepredadores={F_e:.1f}",
xy=(R_e, F_e), xytext=(R_e+80, F_e+7), fontsize=12,
arrowprops=dict(arrowstyle="->", color="black"))
plt.xlabel('Presas')
plt.ylabel('Depredadores')
plt.title('Espacio de fase con método de RK2 (presa v/s depredador)')
plt.subplot(3, 2, 5)
plt.plot(results_rk4['times'], results_rk4['preys'], label='Presas (RK4)', color='blue')
plt.plot(results_rk4['times'], results_rk4['predators'], label='Depredadores (RK4)', color='red')
plt.xlabel('Tiempo')
plt.xlim(0, 2_000)
plt.ylabel('Población')
plt.legend(loc='upper right')
plt.title('Método de Runge-Kutta 4')
plt.subplot(3, 2, 6)
plt.plot(results_rk4['preys'], results_rk4['predators'], color='purple')
plt.axvline(R_e, color='blue', linestyle='--', alpha=0.7, label='Equilibrio Presas')
plt.axhline(y=F_e, color='red', linestyle='--', alpha=0.7, label='Equilibrio Depredadores')
plt.scatter(R_e, F_e, color='black', zorder=5)
plt.annotate(f"Equilibrio:\nPresas={R_e:.1f}\nDepredadores={F_e:.1f}",
xy=(R_e, F_e), xytext=(R_e+80, F_e+7), fontsize=12,
arrowprops=dict(arrowstyle="->", color="black"))
plt.xlabel('Presas')
plt.ylabel('Depredadores')
plt.title('Espacio de fase con método de RK4 (presa v/s depredador)')
plt.tight_layout()
plt.show();
Podemos visualizar que los diagramas de fase resultantes de las condiciones inciales establecidas son figuras ovaladas que parecieran estar centradas en un punto. Este punto es uno especial, puesto que una de las soluciones particulares del sistema de ecuaciones diferenciales.
Para determinar dicho valor, es necesario considerar cuál es la solución cuando las derivadas son cero, es decir, no hay un cambio en las poblaciones cuando cambiamos el valor del tiempo.
$$ \begin{cases} \frac{dp(t)}{dt} &= \alpha p(t) - \beta p(t) P(t) \\ \frac{dP(t)}{dt} &= \beta \epsilon p(t) P(t) - \gamma P(t) \end{cases} $$ $$ \begin{cases} \frac{dp(t)}{dt} &= 0 \implies \alpha p(t) - \beta p(t) P(t) &= 0 \implies (\alpha - \beta P(t))p(t) &= 0 \\ \frac{dP(t)}{dt} &= 0 \implies \beta \epsilon p(t) P(t) - \gamma P(t) &= 0 \implies (\beta \epsilon p(t) - \gamma)P(t) &=0 \end{cases} $$ $$ \therefore \left( p(t) = 0 \, \land \, P(t) = 0 \right) \, \lor \, \left( p(t) = \frac{\gamma}{\beta \epsilon} \, \land \, P(t) = \frac{\alpha}{\beta} \right) $$
De este modo, podemos determinar que existen dos puntos de quilibrio, el primero ocurre cuando las poblaaciones de las presas y los depredadores son cero, lo que nos lleva al punto $(0, 0)$. El segundo es cuando las poblaciones son distintas de cero, lo que nos lleva al punto $ \left( \frac{\gamma}{\beta \epsilon}, \frac{\alpha}{\beta} \right) $.
Podemos continuar con el análisis al utilizar el Jacobiano asociado al sistema de ecuaciones diferenciales: $$ J(x, y) = \begin{pmatrix} \alpha - \beta P(t) & -\beta p(t) \\ \beta \epsilon P(t) & -\gamma + \beta \epsilon p(t) \end{pmatrix} $$
Al evaluar con el punto $(0, 0)$, obtenemos el sisguiente Jacobiano $$ J(0, 0) = \begin{pmatrix} \alpha & 0 \\ 0 & -\gamma \end{pmatrix} $$ De esta matriz se obtienen los valores propios $\lambda_1 = \alpha > 0$ y $\lambda_2 = -\gamma < 0$. Esto implica que el punto de equilibrio $(0,0)$ es un punto silla, ya que en una dirección (asociada al valor propio positivo) las soluciones tienden a crecer y se alejan del equilibrio, mientras que en la otra dirección (asociada al valor propio negativo) las soluciones decrecen y se acercan al equilibrio. Desde el punto de vista físico, este comportamiento refleja que la población nula es inestable: cualquier pequeña cantidad de presas puede crecer en ausencia de depredadores, mientras que los depredadores no pueden sobrevivir sin presas.
Ahora, al evaluar con el punto $ \left( \frac{\gamma}{\beta \epsilon}, \frac{\alpha}{\beta} \right)$, obtenemos el siguiente Jacobiano $$ J\left( \frac{\gamma}{\beta \epsilon}, \frac{\alpha}{\beta} \right) = \begin{pmatrix} 0 & -\frac{\gamma}{\epsilon} \\ \epsilon \alpha & 0 \end{pmatrix} $$ De esta matriz, los valores propios son $\lambda_{1,2} = \pm \sqrt{\alpha \gamma}i$, y como son puramente imaginarios, el punto de equilibrio es un punto centro en donde existen elipses o círculos en torno a este. Esta representación de las curvas de nivel físico porque representan el ciclo de la dinámica de las poblaciones de presas y depredadores, en donde si uno aumenta el otro aumenta pero de forma desfasada, asimismo cuando disminuyen las poblaciones.
En el diagrama de fase para el método de Euler se aprecia un grosor característico que no presentan los métodos de $RK2$ y $RK4$. Esto se debe al error acumulado propio del método de Euler, el cual introduce pequeñas desviaciones en cada paso de integración. Estas desviaciones hacen que la trayectoria no siga una órbita cerrada perfectamente definida, sino que se disperse alrededor de la curva esperada, generando un trazo más ancho.
Por el contrario, los métodos de $RK2$ y $RK4$ muestran órbitas mucho más precisas y ajustadas, con curvas bien definidas en el espacio de fases. Esto refleja su mayor orden de convergencia: mientras Euler acumula errores lineales con el paso de tiempo, $RK2$ los reduce cuadráticamente y $RK4$ de manera aún más eficiente, manteniendo trayectorias cercanas al equilibrio oscilatorio verdadero del sistema.
p_min, p_max = 0, 800
P_min, P_max = 0, 200
nx, ny = 120, 80
p_vals = np.linspace(p_min, p_max, nx)
P_vals = np.linspace(P_min, P_max, ny)
p_grid, P_grid = np.meshgrid(p_vals, P_vals)
dpdt = np.zeros_like(p_grid)
dPdt = np.zeros_like(P_grid)
for i in range(p_grid.shape[0]):
for j in range(p_grid.shape[1]):
dpdt[i, j] = prey_poblation(params['a'], params['b'], p_grid[i, j], P_grid[i, j])
dPdt[i, j] = depr_poblation(params['b'], params['m'], p_grid[i, j], P_grid[i, j], params['eps'])
speed = np.hypot(dpdt, dPdt)
speed_norm = speed / (np.nanmax(speed) + 1e-16)
fig, ax = plt.subplots(figsize=(12,8))
c = ax.pcolormesh(p_grid, P_grid, speed, shading='auto', cmap='viridis', alpha=0.6)
fig.colorbar(c, ax=ax, label='Magnitud del campo')
lw = 0.8 + 2.2 * speed_norm
strm = ax.streamplot(
p_vals, P_vals, dpdt, dPdt,
density=1.2,
color=speed,
linewidth=lw,
cmap='plasma',
arrowsize=1.0,
arrowstyle='->',
minlength=0.1
)
ax.set_xlim(p_min, p_max)
ax.set_ylim(P_min, P_max)
ax.set_xlabel('Presas')
ax.set_ylabel('Depredadores')
ax.set_title('Campo de direcciones (Lotka-Volterra)')
ax.grid(alpha=0.3)
plt.show()
En este diagrama, cada flecha indica la dirección en la que evoluciona el sistema a partir de un estado inicial específico de poblaciones. Tanto el grosor de las flechas como el color del mapa de calor representan la magnitud del campo, es decir, la velocidad de cambio. Las flechas más gruesas y los tonos más amarillos corresponden a variaciones más rápidas de las poblaciones, mientras que las flechas más delgadas y los colores más morados indican cambios más lentos. Además, mientras más lejos se esté del punto de equilibrio, el mapa de calor se torna más amarillo y las flechas se vuelven más gruesas, lo que refleja que la magnitud del campo crece con la distancia al equilibrio.
A medida que nos alejamos del punto de equilibrio de coexistencia, el rango entre los valores mínimos y máximos de las poblaciones aumenta. Esto significa que las oscilaciones en la dinámica poblacional respecto al tiempo $t$ se vuelven más amplias: los picos de crecimiento de presas y depredadores son más altos, mientras que en los descensos las poblaciones alcanzan valores mucho más bajos. En otras palabras, cuanto más distante esté el estado inicial del equilibrio, más marcadas serán las fluctuaciones en las trayectorias del sistema.
y0_cercana = [140, 20]
y0_lejano = [700, 150]
# Configuración de la simulación. Paso de tiempo y número total de pasos.
dt = 0.01
T = 50_000
#Para condiciones iniciales cercanas
results_cercano = simulate(params, y0_cercana, dt, T, rk4_step)
results_lejano = simulate(params, y0_lejano, dt, T, rk4_step)
plt.figure(figsize=(16, 6))
# Gráfico 1: condiciones iniciales cercanas al punto de equilibrio
plt.subplot(1, 2, 1)
plt.plot(results_cercano['times'], results_cercano['preys'], label='Presas (cercano)', color='blue', alpha=0.7)
plt.plot(results_cercano['times'], results_cercano['predators'], label='Depredadores (cercano)', color='red', alpha=0.7)
plt.xlabel('Tiempo')
plt.ylabel('Población')
plt.title('Condiciones iniciales cercanas al punto de equilibrio')
plt.legend(loc='upper right')
plt.xlim(0, 2000)
# Gráfico 2: condiciones iniciales lejanas al punto de equilibrio
plt.subplot(1, 2, 2)
plt.plot(results_lejano['times'], results_lejano['preys'], label='Presas (lejano)', color='blue', alpha=0.7)
plt.plot(results_lejano['times'], results_lejano['predators'], label='Depredadores (lejano)', color='red', alpha=0.7)
plt.xlabel('Tiempo')
plt.ylabel('Población')
plt.title('Condiciones iniciales lejanas al punto de equilibrio')
plt.legend(loc='upper right')
plt.xlim(0, 2000)
plt.tight_layout()
plt.show()
5. Análisis de Parámetros¶
Parámetro $\alpha$¶
params = {
'b': 0.003,
'm': 0.01,
'eps': 0.02
}
a_vals = np.linspace(0.001, 0.1, 20)
dt = 0.01
T = 5_000
a_results = {}
alpha_unicode = '\u03B1'
for a in a_vals:
params['a'] = a
results = simulate(params, y0, dt, T, rk4_step)
R_e = params['m'] / (params['eps'] * params['b'])
F_e = params['a'] / params['b']
a_results[a] = {
'results': results,
'R_e': R_e,
'F_e': F_e
}
a_values = np.random.choice(a_vals, size=5, replace=False)
a_values = np.sort(a_values)
fig, axes = plt.subplots(len(a_values), 2, figsize=(16, 5*len(a_values)))
for i, a in enumerate(a_values):
params['a'] = a
res = simulate(params, y0, dt, T, rk4_step)
t = np.array(res['times'])
prey = np.array(res['preys'])
pred = np.array(res['predators'])
R_e = params['m'] / (params['eps'] * params['b'])
F_e = params['a'] / params['b']
# Serie temporal
x_min, x_max = 0, 2000
mask = (t >= x_min) & (t <= x_max)
y_max = max(prey[mask].max(), pred[mask].max()) * 1.1
axes[i, 0].plot(t, prey, label='Presas', color='blue')
axes[i, 0].plot(t, pred, label='Depredadores', color='red')
axes[i, 0].set_xlabel('Tiempo')
axes[i, 0].set_xlim(x_min, x_max)
axes[i, 0].set_ylabel('Población')
axes[i, 0].set_ylim(0, y_max * 1.1)
axes[i, 0].legend(loc='upper right')
axes[i, 0].set_title(f'Serie temporal ({alpha_unicode}={a:.3f})')
# Espacio de fase
axes[i, 1].plot(prey, pred, color='purple')
axes[i, 1].axvline(R_e, color='blue', linestyle='--', alpha=0.7, label='Equilibrio Presas')
axes[i, 1].axhline(F_e, color='red', linestyle='--', alpha=0.7, label='Equilibrio Depredadores')
axes[i, 1].scatter(R_e, F_e, color='black', zorder=5)
axes[i, 1].annotate(f"Equilibrio:\nPresas={R_e:.1f}\nDepredadores={F_e:.1f}",
xy=(R_e, F_e),
xytext=(R_e-20, F_e+2),
fontsize=10,
arrowprops=dict(arrowstyle="->", color="black"),
bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="black", alpha=0.7))
axes[i, 1].set_xlabel('Presas')
axes[i, 1].set_ylabel('Depredadores')
axes[i, 1].set_title(f'Espacio de fase ({alpha_unicode}={a:.3f})')
axes[i, 1].legend(loc='upper right')
plt.tight_layout()
plt.show();
Al variar $\alpha$ observamos que cambia el punto de equilibrio de la población de depredadores (porque $P^{*} = \frac{\alpha}{\beta}$), y esa traslación del equilibrio altera notablemente los ciclos poblacionales. Con $\alpha$ pequeño, el equilibrio de depredadores queda muy cerca de cero y las oscilaciones tienden a ser extremas: las presas alcanzan picos muy altos luego caen a valores muy bajos, con respuestas retardadas y agudas de los depredadores. Ese comportamiento amplificado aumenta el riesgo de extinción local, sobre todo cuando la tasa intrínseca de crecimiento de las presas es menor que la tasa efectiva de interacción presa-depredador. En resumen $\alpha$ actúa como un regulador de la intensidad dinámica: valores menores empujan al sistema hacia colapsos o extinciones, mientras que valores mayores favorecen ciclos más pronunciados pero también poblaciones de depredadores más elevadas, aunque el ciclo de la dinámica poblacional es más estable y más redondo. Por lo tanto, pequeñas variaciones en $\alpha$ pueden determinar si el sistema se estabiliza, oscila moderadamente o entran en un patrón de "auge y caída".
# convertir
!jupyter nbconvert --to html "Tarea1_Final.ipynb"
# abrir el HTML generado en el navegador por defecto (Windows)
!start "Tarea1_Final.html"